Strzeszczenie

Poniższa praca przedstawia analizę zbioru danych śledzie. Wynikiem końcowym jest wniosek przedstawiający domniemaną przyczynę zmniejszenia długości śledzi. W procesie przetwarzania danych przeanalizowano zbiór wejściowy i wyczyszczono dane. Została także przedstawiona macierz korelacji zmiennych, oraz interaktywny wykres zawartości planktonu w wodzie zależnie od aktualnego miesiąca w roku.

Analiza podstawowa

Inicjalizacja

Załadowanie danych i analiza

data <- read.csv("sledzie.csv")
summary(data)
##        X             length        cfin1              cfin2          
##  Min.   :    0   Min.   :19.0   Length:52582       Length:52582      
##  1st Qu.:13145   1st Qu.:24.0   Class :character   Class :character  
##  Median :26291   Median :25.5   Mode  :character   Mode  :character  
##  Mean   :26291   Mean   :25.3                                        
##  3rd Qu.:39436   3rd Qu.:26.5                                        
##  Max.   :52581   Max.   :32.5                                        
##     chel1              chel2              lcop1              lcop2          
##  Length:52582       Length:52582       Length:52582       Length:52582      
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##       fbar             recr              cumf             totaln       
##  Min.   :0.0680   Min.   : 140515   Min.   :0.06833   Min.   : 144137  
##  1st Qu.:0.2270   1st Qu.: 360061   1st Qu.:0.14809   1st Qu.: 306068  
##  Median :0.3320   Median : 421391   Median :0.23191   Median : 539558  
##  Mean   :0.3304   Mean   : 520367   Mean   :0.22981   Mean   : 514973  
##  3rd Qu.:0.4560   3rd Qu.: 724151   3rd Qu.:0.29803   3rd Qu.: 730351  
##  Max.   :0.8490   Max.   :1565890   Max.   :0.39801   Max.   :1015595  
##      sst                 sal            xmonth            nao          
##  Length:52582       Min.   :35.40   Min.   : 1.000   Min.   :-4.89000  
##  Class :character   1st Qu.:35.51   1st Qu.: 5.000   1st Qu.:-1.89000  
##  Mode  :character   Median :35.51   Median : 8.000   Median : 0.20000  
##                     Mean   :35.51   Mean   : 7.258   Mean   :-0.09236  
##                     3rd Qu.:35.52   3rd Qu.: 9.000   3rd Qu.: 1.63000  
##                     Max.   :35.61   Max.   :12.000   Max.   : 5.08000

Zauważamy że kolumny, które powinny być typu numerycznego są ciągiem znaków.

Próbka danych

head(data)
##   X length   cfin1   cfin2   chel1    chel2   lcop1    lcop2  fbar   recr
## 1 0   23.0 0.02778 0.27785 2.46875        ? 2.54787 26.35881 0.356 482831
## 2 1   22.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 3 2   25.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 4 3   25.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 5 4   24.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 6 5   22.0 0.02778 0.27785 2.46875 21.43548 2.54787        ? 0.356 482831
##        cumf   totaln           sst      sal xmonth nao
## 1 0.3059879 267380.8 14.3069330186 35.51234      7 2.8
## 2 0.3059879 267380.8 14.3069330186 35.51234      7 2.8
## 3 0.3059879 267380.8 14.3069330186 35.51234      7 2.8
## 4 0.3059879 267380.8 14.3069330186 35.51234      7 2.8
## 5 0.3059879 267380.8 14.3069330186 35.51234      7 2.8
## 6 0.3059879 267380.8 14.3069330186 35.51234      7 2.8

Zauważamy że miejscami pojawia się znak ? zamiast danych. W następnym kroku będzimy dążyli do zamiany go na NA, oraz uzupełnienie średnią kolumny.

Czyszczenie danych

#zamień wszystkie kolumny na liczbę rzeczywistą
data2 <- as.data.frame(lapply(data, as.numeric))
#zamień wszystkie wartości NA na średnią dla kolumny
data2 <- data2 %>% mutate_if(is.numeric, ~replace_na(.,mean(., na.rm = TRUE)))
#zamień wszystkie kolumny będące wcześniej liczbą całkowitą na liczby całkowite
data2 <- transform(data2, recr = as.integer(recr), xmonth = as.integer(xmonth))
#wyrzuć kolumnę id
data2 <- subset(data2, select = -c(X))

colSums(is.na(data2))
## length  cfin1  cfin2  chel1  chel2  lcop1  lcop2   fbar   recr   cumf totaln 
##      0      0      0      0      0      0      0      0      0      0      0 
##    sst    sal xmonth    nao 
##      0      0      0      0

Podsumowanie zbioru po czyszczeniu

summary(data2)
##      length         cfin1             cfin2             chel1       
##  Min.   :19.0   Min.   : 0.0000   Min.   : 0.0000   Min.   : 0.000  
##  1st Qu.:24.0   1st Qu.: 0.0000   1st Qu.: 0.2778   1st Qu.: 2.469  
##  Median :25.5   Median : 0.1333   Median : 0.7012   Median : 6.083  
##  Mean   :25.3   Mean   : 0.4458   Mean   : 2.0248   Mean   :10.006  
##  3rd Qu.:26.5   3rd Qu.: 0.3603   3rd Qu.: 1.9973   3rd Qu.:11.500  
##  Max.   :32.5   Max.   :37.6667   Max.   :19.3958   Max.   :75.000  
##      chel2            lcop1              lcop2             fbar       
##  Min.   : 5.238   Min.   :  0.3074   Min.   : 7.849   Min.   :0.0680  
##  1st Qu.:13.589   1st Qu.:  2.5479   1st Qu.:17.808   1st Qu.:0.2270  
##  Median :21.435   Median :  7.1229   Median :25.338   Median :0.3320  
##  Mean   :21.221   Mean   : 12.8108   Mean   :28.419   Mean   :0.3304  
##  3rd Qu.:27.193   3rd Qu.: 21.2315   3rd Qu.:37.232   3rd Qu.:0.4560  
##  Max.   :57.706   Max.   :115.5833   Max.   :68.736   Max.   :0.8490  
##       recr              cumf             totaln             sst       
##  Min.   : 140515   Min.   :0.06833   Min.   : 144137   Min.   :12.77  
##  1st Qu.: 360061   1st Qu.:0.14809   1st Qu.: 306068   1st Qu.:13.63  
##  Median : 421391   Median :0.23191   Median : 539558   Median :13.86  
##  Mean   : 520367   Mean   :0.22981   Mean   : 514973   Mean   :13.87  
##  3rd Qu.: 724151   3rd Qu.:0.29803   3rd Qu.: 730351   3rd Qu.:14.16  
##  Max.   :1565890   Max.   :0.39801   Max.   :1015595   Max.   :14.73  
##       sal            xmonth            nao          
##  Min.   :35.40   Min.   : 1.000   Min.   :-4.89000  
##  1st Qu.:35.51   1st Qu.: 5.000   1st Qu.:-1.89000  
##  Median :35.51   Median : 8.000   Median : 0.20000  
##  Mean   :35.51   Mean   : 7.258   Mean   :-0.09236  
##  3rd Qu.:35.52   3rd Qu.: 9.000   3rd Qu.: 1.63000  
##  Max.   :35.61   Max.   :12.000   Max.   : 5.08000

Korelacje między zmiennymi

corelation <- cor(data2)
corrplot(corelation,method="color", order="hclust", sig.level = 0.01)

Zagęszczenie planktonu w poszczególnych miesiącach

#na podstawie: https://stackoverflow.com/questions/68471610/r-plotly-dropdown-variable-selection-with-color

plot_ly(data = data2, x = ~factor(xmonth), y = ~cfin1, type ="box", xaxis = list(title = "Miesiąc")) %>%
  layout(
    annotations = list(
      list(
        text = "<b>Rodzaj planktonu: </b>", x=0.04, y=1.13, xref='paper', yref='paper',xanchor = "left", showarrow=FALSE
      )
    ),
    xaxis = list(title = 'Miesiąc'),
    updatemenus = list(
      list(
        type = "list",
        x = 0.25,
        xanchor = "left",
        y = 1.15,
        buttons = list(
          list(
            method = "update",
            args = list(list(y = list(data2$cfin1)),
                        list(yaxis = list(title = "cfin1"))
                        ),
            label = "cfin1"
              ),
          list(
            method = "update",
            args = list(list(y =list(data2$cfin2)),
                        list(yaxis = list(title = "cfin2"))
                        ),
            label = "cfin2"
              ),
          
          list(
            method = "update",
            args = list(list(y = list(data2$chel1)),
                        list(yaxis = list(title = "chel1"))
                        ),
            label = "chel1"
              ),
          
          list(
            method = "update",
            args = list(list(y = list(data2$chel2)),
                        list(yaxis = list(title = "chel2"))
                        ),
            label = "chel2"
              ),
          
          list(
            method = "update",
            args = list(list(y = list(data2$lcop1)),
                        list(yaxis = list(title = "lcop1"))
                        ),
            label = "lcop1"
              ),
          
          list(
            method = "update",
            args = list(list(y = list(data2$lcop2)),
                        list(yaxis = list(title = "lcop2"))
                        ),
            label = "lcop2"
              )
        )
      )
    )
  )

Przewidywanie rozmiaru śledzia

Podział zbioru na 2 części

inTraining <- 
    createDataPartition(
        y = data2$length,
        p = .80,
        list = FALSE)

training <- data2[ inTraining,]
testing  <- data2[-inTraining,]

Wybór parametrów algorytmu

metoda repeatedcv (powtórna walidacja krzyżowa) dzieli zbiór treningowy na zbiór treningowy i walidacyjny.

ctrl <- trainControl(
    method = "repeatedcv",
    classProbs = TRUE,
    number = 2,
    repeats = 5)

Utworzenie modelu

fit <- train(length ~ .,
             data = training,
             method = "rf",
             metric = 'RMSE',
             trControl = ctrl,
             ntree = 15)
## Warning in train.default(x, y, weights = w, ...): cannnot compute class
## probabilities for regression
fit
## Random Forest 
## 
## 42067 samples
##    14 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (2 fold, repeated 5 times) 
## Summary of sample sizes: 21034, 21033, 21033, 21034, 21033, 21034, ... 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE      
##    2    1.167790  0.5000386  0.9234103
##    8    1.167366  0.5017007  0.9185563
##   14    1.175046  0.4960310  0.9241144
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 8.

Wybrano drzewo mtry=2, z RMSE = 1.16 i R^2 = 0.501. Są to wartości mocno nieidealne, ale w tym momencie najlepsze.

Predykcja

pred <- predict(fit, testing)
postResample(pred = pred, obs = testing$length)
##      RMSE  Rsquared       MAE 
## 1.1513554 0.5184248 0.9037917

Powód zmniejszania rozmiaru śledzi

Używamy metody varImp do zbadania istotności zmiennych.

imp <- varImp(fit)
plot(imp)

Długość śledzi zmniejsza się głównie przez temperaturę przy powierzchni wody, oraz rozmiary połowów. Najmniej na długość śledzia wpływa miesiąc połowu, ilość pozostawionych śledzi przy połowie i zagęszczenie planktonu.